Simultaneous Estimation of Dimension, States and Measurements: 
Computation of representative density matrices and POVMs 
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This investigation continues a program aiming at obtaining effective quantum models to describe 
measurement statistics. In [Stark, arXiv:1209.5737 (2012)], we have described how the Gram matrix 
G associated to the prepared states and the performed POVM elements can be estimated via convex 
relaxation of a rank minimization problem. This Gram matrix G determines the density matrices and 
the POVM elements uniquely up to simultaneous rotations (with respect to the Hilbert- Schmidt 
inner product) in the vector space of Hermitian matrices. However, when the description of the 
experiment needs to be connected to textbook quantum mechanics, explicit expressions for the 
states and the POVM elements in terms of positive semidefinite matrices are required. In this 
paper, we describe a heuristic algorithm that takes the state-meaurement Gram matrix G as input 
and searches explicit realizations of G in terms of quantum mechanically valid density matrices and 
POVM elements. 



Assume you have built up an experiment that allows 
you to prepare a quantum system in different states, and 
subsequently perform different measurements on this sys- 
tem. The catch is that you do not know exactly what 
states are actually prepared, and what the measurement 
devices truly measure. But maybe you are willing to as- 
sume that the experiment can be repeated such that one 
can gather statistics. Then, the question is what can be 
deduced from this knowledge. We have previously intro- 
duced a general program for simultaneous estimation of 
states and measurements. In [l[ we described how the 
Gram matrix, G, associated to the prepared states and 
the POVM elements can be estimated via convex relax- 
ation of a rank minimization problem. We showed that G 
is not always uniquely determined by the acquired mea- 
surement data. In [2| we introduced a method to test 
whether or not specific assumptions on the Gram matrix 
can guarantee its uniqueness. The Gram matrix G only 
determines the density matrices and the POVM elements 
uniquely up to simultaneous rotations (with respect to 
the Hilbert-Schmidt inner product) in the vector space 
of Hermitian matrices. If we compute them via stan- 
dard factorizations of G, we typically do not get matri- 
ces which are positive semidefinite as we normally would 
expect of proper density operators and POVM elements. 
In this paper we complete the proof of principle of the 
simultaneous estimation of dimension, states and mea- 
surements, in terms of a heuristic algorithm that delivers 
a collection of positive and normalized density operators 
for the states, and positive POVM elements for the mea- 
surements. 

We observe a growing interest in connection with inad- 
equate quantum modeling: In Q the authors introduce 
a method which allow for the falsification of inaccurate 
modelling of the performed POVMs. Temporal drifts 
in the measurement setup are discussed in Q. The pa- 
pers investigate consequences of inappropriate theo- 
retical models. In our work, dimension estimation is tied 
to the estimation of the states and the measurements. 
In 0-[ll| on the other hand, dimension is estimated via 
linear witnesses. 



This paper is organized as follows. In Section U we in- 
troduce the general setting, define the general task more 
precisely, and clarify its role in the general program out- 
lined in [l[ and @ . In Section[H]we discuss the main ideas 
behind our algorithm, and introduce its different subrou- 
tines. Section IlIII considers the necessary 'monitoring' of 
the evolution of the computation, and how this is used for 
the switching between the different subroutines. In Sec- 
tion |IV]wc provide some numerical examples to illustrate 
the application of this technique. 



I. SETTING 

We imagine that we are able to prepare W differ- 
ent states (pw)w and to measure V choices of POVMs 
(E v k)vk (k enumerates the different outcomes). To keep 
the notation simple, we assume that the number of out- 
comes, K, is the same for all measurements. Both the 
density matrices and the POVMs are unknown a pri- 
ori even though you might have an educated guess what 
they might look like. Performing independent repetitions 
of each of the measurements on each of the states, we 
can determine frequency distributions (f WjV k)k for each 
pairing of a state 'u>' with a measurement l v\ Here, 
fw,vk denotes the frequency for measuring 'fc' in case we 
have prepared state 'w' and measured POVM V. In the 
asymptotic limit (i.e., the number of measurement rep- 
etitions goes to infinity), f w . v k is equal to its associated 
probability p^^fc. Consequently, 



fw.vk = te(E v kPw) = 2^ (E v k)ij{Pw)ij — {E vk ) T p w - 

-* ,2 

Here, we have created the vectors E v k and p w in C 
out of the matrices E v k , p w by stacking all the matrix- 
columns on top of each other. When collecting all the 
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frequencies within one data table 

/ A, (1,1) ' ' ' fl,(V,K) 

V = : : 

\ /w,(l,l) ' ' " fw,{V,K) 

and all the vectors p w and E v k within one matrix 
P = (Pi I ■ • • I Pw I E lx | • ■ ■ Bvif), 



(1) 



(2) 



we observe that T> appears as off-diagonal blocks in the 
state-measurement Gram matrix G: 



G = P 1 P 





V 







(3) 



Guided by Occam's razor, we described in how G can 
be estimated via the convex relaxation of a rank mini- 
mization problem. In Q we have studied the extent to 
which the state-measurement Gram matrix is evidence 
based, i.e., to which extent G is determined by the mea- 
surement data T> alone. The state-measurement Gram 
matrix determines the density matrices and the POVM 
elements uniquely up to simultaneous rotations (with re- 
spect to the Hilbert-Schmidt inner product) of all these 
Hcrmitian matrices. 

However, in some applications knowing G might not be 
sufficient and explicit realizations of the states and the 
POVM elements in terms of quantum mechanically valid 
matrices (positive semidefinite) are desired. Finding ex- 
plicit, quantum mechanically valid, representations is the 
goal of this paper. For this purpose, we have designed 
a heuristic algorithm which takes the state-measurement 
Gram matrix G as input, and returns proper states and 
POVMs. From a geometric perspective, our task is to 
fit W + VK (i.e., total number of states plus total num- 
ber of POVM elements) matrices into the cone of positive 
semidefinite matrices such that their pairwise inner prod- 
ucts satisfy the constraints described by G. Thus, the 
wanted states-measurement matrix P can be regarded as 
the solution to the optimization problem 



argmin \\P T P - G\\j 
subject to Columns of P carry valid 
states and measurements. 



(4) 



Here, the constraints refer to the usual demands p w > 0, 
tr{p w ) = 1 for all states, and E v k > 0, J2k Evk = I for 
all POVM elements. Unfortunately, this is not a convex 
optimization problem. (Analogous to min^gij \x 2 — 1|.) 
It might thus seem that (|4]) is intractable. However, in 
contrast to many other non-convex optimization prob- 
lems, we have the advantage that we easily can verify 
when we have reached a global minimum, since the op- 
timal P satisfies P T P = G. This holds true as long as 
G is a Gram matrix that can approximately be recon- 
structed by a quantum model. It should be stressed that 
the algorithm cannot be treated as a 'black box solution' 



to the task of finding proper states and POVMs. We 
cannot exclude that its application may require modi- 
fications of the suggested parameter values in the algo- 
rithm. If the algorithm does not converge quickly enough 
to a desired realization, the evolution of the computation 
needs to be analyzed to adjust these parameters. This 
applies mainly to the parameters dictating the switch- 
ing from one phase of the algorithm (REGULAR, SE- 
LECTION.OF.FASTEST, and PARTIAL; see below) to 
another phase. 



II. THE ALGORITHM 

Here we present the an algorithm to solve (d]). The 
main idea is to iteratively solve sequences of convex op- 
timization problems. To this end we first observe that 
although non-convex, problem (UJ) has a lot of structure. 
When the first factor is replaced by a constant matrix, 
P T i y Pq , we arrive at a family of independent least 



squares problems, since 



\Pn P G\\l = 



j 



\PTP(:,j)-G(:,j)\\1 (5) 



Here, P(:,j) and G(:,j) refer to the j-th columns of P 
and G respectively. Each summand of the RHS of (([5])) 
corresponds to a convex optimization problem, which 
thus can be solved efficiently and reliably. This obser- 
vation leads to Algorithm [TJ see Listing [2] for the descrip- 
tion of the routine 'REGULAR'. There, mat(-) refers to 
the inverse of the transition p i— > p . 

Algorithm 1 Draft of SESAM_rcalizations 
Require: Gram matrix G, Hilbert space dimension d, error 
threshold 

1: Initialize P (cf. Eq. ©) with random density matrices 

and POVM elements. 
2: error := oo 

3: while error > threshold do 
4: for ind = l:size(P,2) do 
5: call REGULAR 

6: end for 

7: error = ||P T P-G|| 2 /||G'||2 
8: end while 

9: Density matrices and POVM elements correspond to 
columns of P; cf. Eq. ([2]). 



Listing 2 REGULAR 
Determine P(:,ind) via 

argmimj \\P T v - G : ,i n d||2 

subject to mat(v) positive semidefinite 



By explicit tests of Algorithm [T] on simple models, one 
observes that it generally appears to work well on single 
qubits, or when we remove the constraints that mat(P(: 
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,j)) must be positive semidefinite for all j. However, 
for estimations in (d > 2)-dimcnsional quantum systems, 
the Algorithm Q] appears to easily get stuck away from 
a global minimum. We suspect that this is connected 
to the somewhat 'edgy' nature of the cone of positive 
matrices over <C d (d > 2). 

The observation that the vectors v in subroutine REG- 
ULAR get stuck reflects that the demands imposed on v 
are too restrictive if mat (if) is cornered in an unfortunate 
part of the cone of positive semidefinite matrices. This 
brings us to the introduction of subroutine PARTIAL 
(see Algorithm [31 line [TU] and Listing [1]). 

Algorithm 3 SESAM_rcalizations 
Require: Gram matrix G, Hilbert space dimension d, error 
threshold 

1: Initialize P (cf. Eq. ([2j) with random density matrices 
and POVM elements, 
mode := 'PARTIAL' 
error := oo 

Use G to compute euclidean distances between columns 
of P. Store all relative distances within a matrix D le \; 

D r el{i,j) ■— y Gu — 2Gij + Gjj 

while error > threshold do 
for ind = l:size(P,2) do 
okLvect := P(:,ind) 
switch mode do 
case 'PARTIAL' 
call PARTIAL 
break 

case 'SELECTION_OF_FASTEST' 
call SELECTION_OF_FASTEST 
break 
case 'REGULAR' 
call REGULAR 
break 
call monitoring 
end for 

call initialization_of_switching_parameters 
call switching_decision 
end while 

Density matrices and POVM elements correspond to 
columns of P; cf. Eq. J2j). 
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Listing 4 PARTIAL 



Use D r ei (cf. line [4] in Algorithm [3} to compute the index 
set ind_NH containing the indices of the d 2 columns of P 
which are closest to P(:,ind). Then, determine P(:,ind) via 

argmin- || [P T ?7]ind_NH - C?ind-NH,ind||2 
subject to mat(7j) positive semidefinite 



The purpose of running PARTIAL is to remove some 
of the demands specified in terms of G(:,j) in the objec- 
tive function in REGULAR. More precisely, PARTIAL 
does not try to align the dynamic P-column v as well 
as possible with respect to all states and measurements. 
Instead it only tries to align v as well as possible with 



the vectors of its local neighborhood. Thus, PARTIAL 
aims at restoring local neighborhoods while REGULAR 
tries to directly arrange the whole family of states and 
measurements. In Listing 0] we chose the minimal choice 
d 2 for the size of the local neighborhood. If we selected 
fewer vectors, the result of each execution of PARTIAL 
would not be unique, since the vectors P(:, j) live in a d 2 - 
dimensional space (recall that the matrices mat(P(:,j)) 
are forced to be Hermitian). 

To run PARTIAL, we need to determine the local 
neighborhood of each column P(:,j) by selecting the d 2 
nearest neighbors (with respect to Hilbert-Schmidt dis- 
tance) of P(:, j). The necessary knowledge about all the 
pairwisc relative distances -D r ci(i,j) between columns i 
and j of P can be computed directly via the Gram ma- 
trix G: 



D lel (i,j) 2 = (P(:,i)-P(:,j),P(:,i) - P(:,j)) 



G 



(6) 



13 • 



This calculation is performed on line [4] of Algo- 
rithm [3] Looking at Algorithm [3] we recognize 
that — for the moment we ignore the routine SELEC- 
TION_OF_FASTEST — in some iterations of the while- 
loop we run REGULAR and in some other iterations we 
run PARTIAL. Thus, we need to prescribe which routine 
to select. This is specified in Listing [5] but we postpone 
the explanation of the switching criteria for later. 

In summary we have arrived at an Algorithm that 
switches back and forth between REGULAR and PAR- 
TIAL. PARTIAL tries to restore local neighborhoods 
while REGULAR tries to find the global, complete con- 
figuration. REGULAR succeeds if it starts with an ini- 
tial condition that is close to a solution. This might be 
achieved via an educated guess. However, even if PAR- 
TIAL manages to arrange all the vectors such that sup- 
posedly near vectors are near, we have not necessarily 
constructed a good starting point for REGULAR. This 
is due to the complicated geometry of the cone of positive 
semidefinite matrices in case of Hilbert space dimensions 
d > 2. The edged shape of the cone of positive semidef- 
inite matrices makes it necessary to rotate the complete 
collection of states-measurement vectors as a whole so 
that the edges of the contour of the states-measurement 
configuration fits into the edges of the cone of positive 
semidefinite matrices. Since we only allow one vector to 
move in each iteration (that is the precondition for using 
convex optimization in each step) this becomes a tiresome 
business because the dynamic vector is held back by all 
the static vectors due to the interactions described by 
G(:,j). The dynamics of the whole states-measurement 
configuration thus resembles the movement of a worm: 
each segment of the worm corresponds to one column 
P(:,j). We need the worm to travel from A to B. Thus, 
each segment of the worm has to travel from A to B. In 
each time step, only one segment of the worm is dynamic 
and wants to reach B but its movement is highly con- 
strained by its bonding to the static segments. Although 
this analogy is admittedly a little far fetched, it neverthe- 
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less provides an intuition for the dynamics of the bunch 
of state-measurement vectors. Speeding up the motion 
of the states-measurements worm is the main purpose of 
the remaining routine. 

The routine SELECTION_OF_FASTEST, see List- 
ing [51 takes into account that some worm segments 
move faster than others, and it operates by ignoring 
the slowest moving segments. Thus, instead of letting 
the worm being slowed down by its small segments, we 
ignore the slow segments. Consequently, during SE- 
LECTION.OF.FASTEST, the evolution of the columns 
P(:,j) only tries to align the dynamic vector with respect 
to the fast moving components of the state-measurement 
bunch. 



Listing 5 SELECTION_OF_FASTEST 
Use shifts (cf. line [2] in Listing [6]) to determine the index 
set ind jast containing the indices of the d 2 columns of P 
which have moved the fastest in the preceding round. Then, 
determine P(:,ind) via 

argmin ff || [P T i7]ind_fa S t - GindJasWndlh 
subject to mat(v) positive semidefinite 



III. SWITCHING SUBROUTINES 

In the previous section we described the main parts 
of the algorithm. Here we describe the criteria for the 
switching between these subroutines. This is specified in 
Listing [U To detect the right moment to switch from 
one routine to the other, we require some basic mon- 
itoring of the computation. The monitoring in List- 
ing [6] records a history of how fast the individual vec- 
tors move (— > shifts), a history of the temporary error 
\\P T P — G|| 2 /||G||2 (— ► error Jiistory) , and the current 
trend of the error (— > trend). Moreover, we store the 
quantities zenith and counter_no_zenith. The variable 
zenith stores the largest error that has been encountered 
until the last switching. The quantity counter_no_zcnith 
counts for how many iterations the procedure has not 
updated the value of zenith. 

Listing [7] initializes some parameters necessary for 
switching_decision. We have chosen random initializa- 
tions to break periodicities (we observed that in some 
situations the complete states-measurement vectors un- 
dergoes periodic movements). 

We have chosen to switch between the different rou- 
tines according to the following sequential ordering: 

REGULAR Q PARTIAL 

(D\ I© 

SELECTION_OF_FASTEST 

(T) We decide to leave REGULAR whenever staying 
in this phase does not lead to meaningful improvements. 



Listing 6 monitoring 

1: shift = ||old_vect - P(:,ind)|j 2 /||P(:,ind)||2 

2: shifts = [shifts shift] 

3: error = ||P T P - G|[ 2 /||G||2 

4: errorjiistory = [raw_history error] 

5: trend = (error - error_history(end-size(P, 2) + l) / er- 

ror_history(end-size(P, 2+1) 
6: if error > zenith then 
7: zenith = error 
8: counter_for_no_zenith = 
9: else 

10: counter_no_zenith = counter_no_zenith + 1 
11: end if 



The variable trend is updated in Listing[6] The threshold 
switch_REG_to_PARTIAL is randomly updated in List- 
ing E 

(2) This transition is enforced under two circum- 
stances. In the first case, we switch whenever 
||shifts(end-2*size(P,2)+l:end)[| 00 < 0.002, i.e., we 
switch whenever PARTIAL gets stuck. The the second 
case occurs when the procedure has not encountered a 
new zenith in no_zenith_threshold-many iterations. This 
rule is motivated as follows: Recall that we start PAR- 
TIAL after REGULAR got stuck in a local minimum 
(eventually just a plateau). We apply PARTIAL with 
the intention that it should lift the states-measurement 
configuration out of the local minimum. To do so, we 
have to cross an 'error-barrier'. If the procedure has 
not reached a new zenith in no_zenith_threshold-many 
iterations, then we take this as an indication that the 
states-measurement configuration may have crossed a 
barrier, and we could check whether or not the bottom 
of the current valley coincides with the global minimum 
||P T P-G|| 2 = 0. 

(3) Before exploring the valley with REGULAR, we 
first apply SELECTION.OF.FASTEST, since we do an- 
ticipate the 'worm-like' dynamics described in the pre- 
vious section. This phase is left when it yields no 
more improvements. Usually, this happens rather quickly 
because typically, SELECTION.OF.FASTEST leads to 
chaotic movements soon after its initiation. 



IV. EXAMPLES 

The goal of this section is to present our experiences 
when running the proposed algorithm [3] in different sit- 
uations. We are considering 2- to 4-dimensional systems 
with different levels of randomness in the states and the 
POVM elements. We cover the following three cases: 

• Scenario 'pure'. We are sampling pure states and 
projective measurements uniformly from the Haar 
measure. 

• Scenario 'partly mixed'. First, we uniformly sample 
r] w from [771b, ?7ub] for each state l w\ and (x v k from 
[/•tibj Mub] for each POVM element i (vk)\ We chose 
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Listing 7 initialization_of_switching_paramcters 

1: 
2: 
3: 
4: 
5: 
6: 
7: 



and 



9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 



Throw a random coin; prob('heads') = 0.7 
if coin='heads' then 

switch_FAST_to_REG = -0.08 
else 

switch_FAST_to_REG = -0.05 
end if 

Throw a random coin; prob('heads') = 0.7 
if coin='heads' then 
9: switch_REG_to_PARTIAL = -0.01 
10: else 

11: switch_REG_to_PARTIAL = -0.02 
12: end if 

13: Throw a random coin; prob('heads') = 0.7 
14: if coin='heads' then 
15: no_zenith_threshold = 3 * size(P, 2) 
16: else 

17: no_zenith_threshold = 7 * size(P, 2) 
18: end if 



Listing 8 switching_decision 

1: 
2: 
3: 



switch mode do 
case 'PARTIAL' 

if ||shifts(end-2*size(P, 2) + l:end)||, 



< 0.002 



then 



4: 
5: 
6: 
7: 



mode = 'SELECTION.OFJFASTEST' 
zenith = -1 
break 
end if 

if counter jor_no_zenith > no_zenith_threshold 



then 



mode = 'SELECTION_OF_FASTEST' 
zenith = -1 
break 
end if 

case 'SELECTION_OF_FASTEST' 

if trend > switch_FAST_to_REG then 

mode = 'REGULAR' 

zenith = -1 

break 
end if 
case 'REGULAR' 

if trend > switch_REG_to_PARTIAL then 

mode = 'PARTIAL' 

zenith = -1 

break 
end if 



Vib = Pib = 0.6 and 7y u b = p u b = 0.8. To generate 
the states, we sample separately for each state a 
unitary U w from the Haar measure, and define 

p w := U w di&g(r] w ,l - r) w , 0, 0) U* . 

To generate the POVM elements, we set P vk , such 
that 



Then, for k = 1,..., (K - 1), 

E v k : = U v k P v k U v f c , 



K-l 



E, 



K 



/c=i 



Here, all the unitary matrices are again sampled 
according to the Haar measure. 

• Scenario 'purified'. All POVM effects on H asso- 
ciated to the K possible outcomes are of the form 
{(k\U\a)}^ =1 where U is a unitary on H ® "Hanc, 
% anc a iVdimensional ancilla system, {\k)}^ =1 is 
an orthonormal basis in 'Hanc, and \a) is an arbi- 
trary state in % anc . Moreover, for all choices of U, 
the definitions 



M, 



\U\a) 



M£M k de- 



yield valid effects and therefore, 
fines a POVM. In the present scenario 'purified', 
we construct the POVM V by sampling U v on 
H ® Hanc with respect to the Haar measure, by 
defining M v k := (k\U v \l), and by setting E v k := 
M* k M v k- To sample the states, we set 

:= diag(p w i, ...,p wd ) 
with p wa sampled uniformly from [0, 1]. Then, 
77 fS ') rj* 

V w Pw u w 



with U w sampled from the Haar measure. 

We have run each of the scenarios 'pure', 'partly 
mixed', and 'purified' 100 times for 2-, 3-, and 4- 
dimensional systems. The stopping criteria was 



max{|tr( piEj) - G tJ \} < 10" 



(7) 



Each run converged successfully. However, we ob- 
serve that some cases required considerably more it- 
erations (i.e., calls of the subroutines PARTIAL, SE- 
LECTION.OF.FASTEST, and REGULAR) than others. 
This is visible in Fig. [1] to Fig. [9l Table U summarizes 
our findings. Calculations have been performed using 
CVX Q1II3 calling SeDuMi 0. 



V. CONCLUSION 

We have introduced a heuristic algorithm to generate 
proper density operators and POVMs from the knowl- 
edge of the Gram matrix of a collection of a priori un- 
known states and measurements. In other words, based 
only on the knowledge of the (Hilbert-Schmidt) inner 
products between states and POVM elements, between 
states and states, and between POVM elements and 
POVM elements, this technique delivers proper density 
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TABLE I: Numerical findings 



Scenario 


d 


W 


V 


A" 


successes 


failures 


'pure' 


2 


8 


4 


2 


100 





'partly mixed' 


2 


8 


4 


2 


100 





'purified' 


2 


8 


4 


2 


100 





'pure' 


3 


27 


9 


3 


100 





'partly mixed' 


3 


27 


9 


3 


100 





'purified' 


3 


27 


9 


3 


100 





'pure' 


4 


56 


9 


4 


100 





'partly mixed' 


4 


56 


9 


4 


100 





'purified' 


4 


56 


9 


4 


100 






estimate the Gram matrix G associated to the involved 
states and POVM elements. In Q we focused on the 
freedom in choosing G, and introduced a method to test 
whether or not specific assumptions on the Gram matrix 
can guarantee its uniqueness. By the extraction of states 
and POVMs presented here, we have thus completed the 
task of simultaneous estimation of dimension, states and 
POVMs. 

However, several questions remain concerning all the 
steps of the estimation procedure, ranging from scenarios 
involving non-asymptotic sampling (leading to statistical 
fluctuations in the Gram matrix G) to efficient, stable, 
and practically applicable numerical procedures. 



operators and POVM elements, i.e., operators that are 
positive semidefinite and satisfy the standard normaliza- 
tion conditions. We do not know to what extent this 
computation of operators in the positive cone, consistent 
with the acquired data, is a computationally hard prob- 
lem. 

This paper concludes a proof of principle that it is pos- 
sible to simultaneously estimate the dimension, states, 
and measurements in the asymptotic regime. This pro- 
gram was introduced in where we showed how to 
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FIG. 1: Scenario 'pure', (d,W,V,K) = (2,8,4,2). 
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FIG. 2: Scenario 'partly mixed', (d,W,V,K) = (2,8,4,2). 
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FIG. 3: Scenario 'purified', (d,W,V,K) = (2,8,4,2). 
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FIG. 4: Scenario 'pure', (d,W,V,K) = (3,27,9,3). 
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FIG. 5: Scenario 'partly mixed', (d,W,V,K) = (3,27,9,3). 
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FIG. 6: Scenario 'purified', (d,W,V,K) = (3,27,9,3). 
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FIG. 7: Scenario 'pure', (d,W,V,K) = (4,56,9,4). 
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FIG. 8: Scenario 'partly mixed', (d,W,V,K) = (4,56,9,4). 
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FIG. 9: Scenario 'purified', (d,W,V,K) = (4,56,9,4). 



